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Figure 1. Free-energy versus Tk/J — T c k/J, where the Monte Carlo estimate of T c k/J = 1.8163 
has been used. The inset figure shows a blow-up of the intersection of the high- and low-temperature 
curves. 
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Table 3. Details of high-temperature expansion in terms of weak lattice constants. Body of table 
gives lattice constant for 8 to 14 bonds on SC lattice. (The contribution from lower-order graphs is 
given in the Appendix). Column 1 identifies graph type, in the notation of Sykes et al 1966. (a. 6) 
denotes a type a and type b graph with a common vertex, [a, 6] is used to denote types a and b 
as separate components. The '8' denotes (p.p), the 'figure 8' type. Column 2 is the weight, with 
a common factor of q — 1 removed. Here S denotes (q — 1), Q denotes (q — 2), T denotes (q — 3), 
W F = Q{q 2 - 5q + 7) and W 117 = Q{q 3 - 9q 2 + 29q - 32), W 41 = Q{q 3 - 6q 2 + 149 - 13), W 69 = 
QT{q 2 - 4q + 5), W 70 = Q 2 {q 2 - 5q + 8), W 78 = Q{q 3 - 7q 2 + 18q - 17), W 79 = Q 2 {q 2 - 5q + 8), W 93 = 
{QT) 2 ,W m = QT(q 2 — 5g + 7), W M = Q{q 2 - 3g + 3), W = {q 3 - 5q 2 + lOq - 7), W Q =Q{q 2 -2q + 2). 



N \N - 1/N] 

16 0.57804 (0.1990)* 

17 0.57908 (0.2090) 

18 0.57856 (0.2041) 

19 0.57857 (0.2042)* 

20 0.57848 (0.2033) 



0.57822 (0.2008)* 

0.57852 (0.2037)* 

0.57838 (0.2025)* 

0.57855 (0.2040)* 

0.57847 (0.2031) 



\N + 1/N] 

0.57870 (0.2054) 

0.57855 (0.2040) 

0.57856 (0.2041)* 

0.57934 (0.2031)* 



Table 4. q = 3, Dlog Pade approximants to the simple-cubic lattice spontaneous magnetisation, 
giving the location and residue of the "pseudo critical point". 
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Table 2. Coefficients in the high-temperature expansion for defined by equation (1). 



16 



Tables and table captions 



n 




m n 







1 


1 





6 


2 


-3 


2 


10 


6 


-18 


24 


11 


6 


-18 


24 


12 


-12 


42 


-56 


14 


30 


-135 


270 


15 


60 


-270 


540 


16 


-96 


477 


-930 


17 


-132 


648 


-1296 


18 


346 


-1980 


4768 


19 


498 


-2988 


7968 


20 


-636 


4140 


-10560 


21 


-2210 


14052 


36922 


22 


3000 


-21690 


64812 


23 


7344 


-52920 


163440 


24 


-7110 


55020 


-165464 


25 


-25836 


201852 


-659088 


26 


17802 


-162774 


600024 


27 


107450 


-914538 


3278256 


28 


-59358 


555750 


-1980408 


29 


-353376 


3229524 


-12285816 


30 


105944 


-1188327 


5005014 


31 


1342914 


-13301370 


55200864 


32 


-77154 


1402686 


-6062712 


33 


-4995004 


52334268 


-227203096 


34 


-226914 


95751 


1954650 


35 


17383710 


-195398208 


914339736 


37 


-64127562 


761838084 


-3742275288 


38 


-32638848 


359664885 


-1761828642 


39 


231546628 


-2910516786 


15132717432 


40 


160963416 


-1946958399 


10380877350 


41 


-805061298 


10681132140 


-58385376120 


42 


-795051840 


10207745148 


-56515869708 


43 


2914349712 


-40522674258 


232316142012 



Table 1. Coefficients in low-temperature expansions for Ao, M and x, denned by equations (3), (5) 
and (6). 
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An interesting additional check is that our general-g expansion enables us to calculate 
K{jp) the mean number of clusters in the low-density regime of the bond percolation 
problem. Bond percolation can be regarded as the q — > 1 limit of the Potts model 
(Fortuin and Kasteleyn, 1972, or, more accesibly, Wu, 1978). Sykes (1986) give series 
for site clusters which agree with our series, after taking into account the cluster of 
1 site and zero bonds that is included in the Potts model limit, but excluded from 
conventional enumerations of bond percolation. 

The agreement between the series of Sykes and the results derived from the weak 
graph expansion of the Potts model give a useful test of our tabulation because the 
original derivation of the percolation series used a quite different type of expansion. 
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We have also shown that series methods can provide not only qualitative 
information about the nature of a phase transition, but also quantitative estimates of 
critical parameters for first-order, as well as second-order phase transitions. 

Our methods have still provided by far the longest extant high-temperature series, 
and this is needed to determine the nature of the phase transition, as well as to identify 
its location. 

Appendix 1. Comparison with weak-graph expansions 

As noted in section 2, our high-temperature series shows significant differences 
from that given by Straley (1974). Since the program used to produce our series (and 
those of I) is our first use of a site representation for calculating high-temperature 
series, these differences caused us some concern. In order to check our program we 
recalculated the free energy using a conventional weak-graph expansion. 

The graphs that are required are those with no vertices of order 1 and in which the 
removal of one line does not increase the number of components. Each graph requires 
a ^-dependent weighting factor. This can be assigned by assigning each directed bond 
an integer flow in the range 1 to q — 1. At each vertex the sum of inwards flows minus 
outwards flows must be modulo q (see for example, Straley, 1974). The weight is 
the number of ways such flows can be assigned. An alternative algebraic formulation 
is given by Domb (1974). For planar graphs, the weighting is equivalent to q~ x times 
the number of §- colourings of the graph and its exterior (see for example Wu, 1982). 

Table 3 gives details of the expansion. To save space, the lowest order terms are not 
shown in the Table. They are, at fourth order, 3 polygons, at sixth order, 22 polygons 
and at seventh order, 18 theta graphs. These then give the following low-order terms, 
with subsequent terms being given in Table 3: 

± In A(<jr = 1) = 3a; 4 + 22a; 6 - 18a; 7 + . . . , 

In A(2) = 3a; 4 + 22a; 6 + Oa; 7 + . . . , 

In A(3) = 6a; 4 + 44a; 6 + 36a; 7 + . . . , 

In A(4) = 9a; 4 66a; 6 + 108a; 7 + . . . , 

In A(5) = 12a; 4 + 88a; 6 + 216a; 7 .... 

The results of the finite lattice calculations agree with the weak graph expansion 
for In A. Further the q = 2 series (see I) reproduces the known Ising model series. 
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terms, as has been done by Vohwinkel, gives AM = 0.463 at the Monte Carlo value 
of T c , a drop of 7%. A similar decrease is found using the extended series at the other 
value of T c used. 

Monte Carlo estimates of this quantity have been obtained by Gavai et al (1989), 
who find AM = 0.395 ± 0.005, only a little lower than our estimate. 

In Table 4 we show the location and residues of Dlog Pade approximants to the 
magnetization series. They indicate a critical temperature of e J l kTc 0.5785, some 
0.3% above the value we believe to be correct. A "pseudo magnetization exponent" 
around 0.2 is also suggested. As shown in II, this is entirely consistent with the 
existence of a first-order transition. 

3.3. Susceptibility 

As for the specific heat, one expects the susceptibility at a first-order transition to 
be undefined, while having well-defined left- and right-hand limits. As we only have 
low-temperature susceptibility series, we can only find one limit. We find a large but 
finite value of the susceptibility at T~ , notably 13 ± 3. on the low-temperature side. 

4. Discussion of results 

Hamer et al (1992) have studied the quantum Hamiltonian version of the 3-state 
Potts model in (2 + 1) dimensions. They find AM = 0.42 ± 0.02, and a latent heat 
jump of 0.21 or 0.24, depending on the lattice. Thus we see that the three distinct 
methods of study, series analysis and Monte Carlo of the Potts model, and series 
analysis of the quantum analogue of the Potts model, give consistent results. That is, 
they all find a fairly weak first-order transition, with a small latent heat, but quite a 
large jump in the magnetisation. 

We have shown how the finite-lattice method can provide competitive series in 
three-dimensions, though the computational complexity ensures that state-of-the- 
art direct methods, well-programmed, will eventually be superior. And as the 
dimensionality increases, the method becomes steadily worse (Guttmann and Enting, 
1993b). 
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(Alves et al 1991), 1.8164 ± 0.0001 (Fukugita et al 1990), 1.8161 ± 0.0001 (Gavai et al 
1989) and 1.81624 ± 0.00006 (Wilson and Vause, 1987). 

It is clear that these estimates do not all agree within the stated precision. However 
it is also clear that the Monte Carlo estimates are all slightly lower than our estimate. 
The average of the Monte Carlo estimates is taken to be 1.8163, and all subsequent 
analysis will be performed using both our central estimate of T c and the Monte Carlo 
average. Note too that the MC estimate is made under the assumption of a first-order 
transition, while our series analysis makes no such assumption. 

By differentiation, we can readily construct series for the high- and low- 
temperature internal energy. Integration of these gives a latent heat of 0.264 ± 0.011, 
where the error is one standard deviation in the average of the approximants. This 
error swamps the error induced by the uncertainty of the critical temperature. Our 
normalisation of the Hamiltonian has, as a consequence, that the internal energy varies 
between and 2. Other workers use a different normalization, in which the internal 
energy varies between and 1. Therefore their latent heat estimates must be doubled 
to be compared with that given here. The specific heat at T~ is found to be 30 ± 4, 
while T c + was considerably lower, at 11.1 ±0.6. For a first order transition, the specific 
heat is undefined at T c , though the left- and right-hand limits are of course defined. 
By use of Monte Carlo methods , Gavai et al (1989) find a latent heat of 0.16 ± 0.008, 
Alves et al find 0.1606 ± 0.0006 while Gavai and Karsch (1992) find a latent heat of 
0.160 ± 0.07. Our result is some 60% higher than these estimates. More precisely, 
we find E(T~) = 1.151 ± 0.009 and E(T+) = 1.414 ± 0.004. We have no explanation 
for the difference between our results and the Monte Carlo results. However, we note 
that our methods did give the correct latent heat in the two-dimensional case. That 
is, all the exact values lay within the range defined by the set of approximants with 
extreme outliers removed. 

3.2. Magnetisation 

As T~ is approached, the gradient decreases rapidly. Consequently, the 
magnetisation gap depends critically on the estimate of the critical temperature. Using 
the MC value of T c , we find AM = 0.498, while using our estimate of T c , we find 
AM = 0.505. Another factor is the length of the series. Extending the series by 13 
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series terms. In principle any order of differential equations can be used, but first- 
order (m = 1) was mostly used in the current work. Finding the coefficients of Qk 
and P reduces to the solution of a system of linear equations, but this system is often 
ill-conditioned, so that care must be taken in its solution. 

These differential equations were then integrated numerically to obtain estimates 
of the desired physical quantities. In all cases a number (up to 10) of DAs were 
integerated and the results averaged to obtain the means and standard deviations 
shown in the tables and graphs below. All calculations were performed in quadruple 
precision (approximately 34 decimal places), so that all series terms could be 
represented without loss of precision. 

We performed the numerical integration with an extrapolation method of the 
Bulirsch-Stoer type, as described by Hairer (1987, Section II. 9). 

The analysis in II involved applying these methods to square-lattice series for 
comparison with the exact results of Baxter (1973, 1982) to asses the extent to which 
the methods could distinguish between first-order and continuous transitions. We 
showed that we could clearly distinguish the order of the transition in all known cases. 

Simple sequence transformations were used to generate the most appropriate series, 
and hence DA, for numerical integration. In general, if a quantity is believed to behave 
like T h at the origin of integration, it is often useful to transform the series X)i=o 
say, to ^2iLo a iT z /T k so that the transformed function approaches a constant at the 
origin. Thus, for example, we worked with the series for xl z& ■> rather than % itself. 
Similarly, instead of the mag netisation M = 1 - 3z e - lSz 10 - lSz 11 .., we worked with 
M — 1 + Sz e . We now discuss our numerical results in greater detail. 

3.1. Internal energy 

We integrated the internal energy series U(z) from T = and T = oo until they 
crossed at T c . The results are shown in Fig. 1, with the intersection region shown 
as an inset. It is perfectly clear already, from this graph alone, that the transition is 
first-order. This follows from the fact that the gradient is clearly discontinuous. The 
intersection of the two curves gives the critical temperature. From a range of several 
approximants, we find the intersection at kT / J = 1.8168 ± 0.0012, which compares 
with Monte Carlo estimates of 1.8166 ± 0.0002 (Yamagata 1993), 1.816454 ± 0.000032 
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The low-temperature series also disagree with those of Miyashita et al (1979) at order 
z 30 . We have not been able to obtain the full field-dependent corrections to their series, 
but note that their coefficient for z 30 fi 5 should be (q — l) 5 times the corresponding 
Ising coefficient (see Sykes et al , 1965), i.e. it should be 44998| and not 45008 as 
published. As our coefficients agree with those of Bhanot et al and Vohwinkel, we are 
confident that they are correct. 

The most serious discrepancies are between our new series and the high- 
temperature series published by Straley (1974), disagreeing at order v 8 and v 10 . In 
view of the gross disagreements at quite low order, we present, in an appendix to this 
paper, an independent re- calculation of these series for general q, using a conventional 
weak-graph expansion. This expansion confirms our finite lattice calculations and 
also reproduces the series for the mean number of clusters in bond percolation. The 
expansion also agrees with Ising model series from I, but this is a weaker test of 
either our weak-graph expansion or the finite lattice method expansion because many 
graph types have zero contribution for q = 2 and, in particular, only even powers of 
v occur. Nevertheless, this appendix provides useful series for general q, and is likely 
to be useful for other workers as a check on any long series that may subsequently be 
obtained. 

3. Analysis of series 

The series generated as described above were all analysed by the method of 
differential approximants (DA) (Guttmann, 1989, page 83ff). This method generalizes 
Pade approximants by fitting an ordinary differential equation of the form 

m 

1=0 

(with d 1 denoting (^) l ) to the available series terms. Here Qk(x) = X^i^fo Qki% z and 
P(x) = X^°o Vi x% are polynomials. We chose q m o = 1, so that the origin is not a 
regular singular point. This allows integration of the differential equation starting at 
x = 0. (This then corresponds to logarithmic derivative Pade approximants when 
m = 1). For magnetisation series, homogeneous DAs (P = 0) are often most useful. 
Generally, the degrees of Qk and P are chosen to use all (or most) of the available 
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Note that for q > 3 an additional 'transverse' susceptibility can be defined (Straley 
and Fisher, 1973). 

Previously, series expansions for the Potts model on the simple cubic lattice had 
been obtained by Straley (1974) ( for low-temperature In Z, M and % to z 2i and high- 
temperature series quoted to v 10 ) and by Miyashita et al (1979) (In Z with full field- 
dependence to order z 33 and full temperature dependence to fi 11 ) In 1993 Bhanot et al 
extended the free-energy and magnetisation series to 39 terms, and the susceptibility 
series to 35 terms. Vohwinkel (1993) has given the magnetisation series to 56 terms, 
and has (unpublished) similar length series for other thermodynamic properties. 

The description of the finite lattice calculations in I was couched in terms of the 
general g-state Potts model and it was that formalism that was used in the present 
study. As noted in I (see also II, equation 17a and b), the amount of memory required 
increases with q. The present calculations used cuboids with cross-sections of up 
to 3 X 4 sites, giving the low-temperature series correct to z 41 . As described in I, 
comparing the 3x4 approximation for q = 2 to higher-order q = 2 calculations allows 
us to determine the correction required to the 3x4 approximation for all other q 
values. This has enabled us to extend the low-temperature series to z i3 . 

For the high-temperature series, we ran on a 4 X 4 lattice, which required storage of 
900 MB. This gave series correct to v 21 . All programs were run on an IBM 3090/400J 
with 1/2 GB of memory and 2 GB of backing storage. One run was also performed 
on a Cray EL. The runs for the high-temperature series took about 50 hours, the 
low-temperature 3x4 lattice runs substantially less. To extend the low-temperature 
series would currently take 4 GB of memory. However, for low-temperature series, we 
believe that the shadow-lattice method is computationally superior, so there seems 
little incentive to develop this method further for low-temperature Potts series. We 
discuss this point further in Guttmann and Enting (1993b). 

The coefficients for \ n , m n , c n and a n for q = 3 are listed in Table 1. These series 
disagree in several places with previously published series. The low-temperature series 
disagree at z 2i with those published by Straley (1974). This error had previously been 
detected using series obtained by the method of partial generating functions (Enting, 
unpublished) and in fact the error was a typographical mistake; the analysis by Straley 
used the correct series (Straley, personal communication). Apparently the term quoted 
in the appendix to Straley (1974) as 486(y 1 6 +y$)x 2i should have been 496(y? +y%)x 2i . 

6 



on q possible values (denoted '0' to q — 1). An energy AE is associated with each pair 
of interacting sites that are in different spin states, and an energy of applies to pairs 
of interacting sites in the same state. We consider the simple cubic lattice, with each 
site interacting only with its 6 nearest neighbours. Each site not in state '0' has an 
additional field energy H . 

The thermodynamic quantities can be derived from the partition function. We 
choose the normalisation such that the state with all sites in state '0' has zero energy. 
In this normalisation, the partition function is commonly denoted A. 

We work in terms of the expansion variables z = exp(— AE/kT), \l = exp(—H/kT) 
and the high-temperature variable v = (1 — z)/(l + (q — l)z) 

For the simple cubic lattice, the high-temperature expansion for the partition 
function takes the form 

A = q- 2 (l + (q-l)z) 3 Hv), (1) 

with 

= a ^ n = i + 3(g - iy + .... 

n 

For the low-temperature expansion, we use a modified field variable x = 1 — \l and 
truncate at order x 2 so that the partition function is expressed as 

A = A + xA 1 + x 2 A 2 + (2) 

The expansion of the zero-field partition function is written as 

A = Y / *nZ n . (3) 

n 

The internal energy is given by 

U = z^/A„. (4) 

az 

the order parameter by 

M = l + -±-!± = Y i m n z n , (5) 
q - 1 A ^ 

and the susceptibility by 



(Knak Jensen et al , 1979, Nienhuis et al 1981, Blote and Swendsen, 1979, Kim and 
Joseph, 1975, Kogut and Sinclair, 1981, Enting, 1974, Herrman, 1979) claimed to see 
evidence of a first-order transition for the three-dimensional 3-state model, but these 
were based on Monte Carlo analyses on small lattices, or series work on short series. 
In 1988, results of large scale QCD calculations by Bacilieri et al (1988) and by Brown 
and Christ (1988) gave conflicting results, with the former claiming to find evidence 
of a second-order phase transition, while the latter claimed to find evidence of a first- 
order transition. Further Monte Carlo studies by Fukugita et al (1990) supported the 
first-order results, but Cabasino et al (1989) argued that the evidence was equally 
good for a first-order or continuous transition. Brown (1989) argued for a first-order 
transition, a result supported by Gavai and Karsch (1992), who considered the effect 
of additional next-nearest neighbour couplings. Further, Alves et al 1991 also used 
high precision Monte Carlo methods to conclude that the transition was first-order. 
Very recently, Bhanot et al (1993) used a method similar to ours to extend the Potts 
model series, and analysed the series to provide additional evidence for a first-order 
transition. After this work was completed, Vohwinkel (1993) showed how the shadow- 
graph method of Sykes (1965) could be used to extend the series even further than 
we have. However Vohwinkel has only presented the magnetisation series, and no 
analysis. 

In order to distinguish between first-order and continuous transitions we have 
analysed the series using differential approximants (Guttmann, 1989) to integrate the 
series. The studies of the square lattice system in II gave a procedure for determining 
the order of the transition in this way. We have used these ideas to locate the critical 
temperature, and to identify the nature of the transition. 

The layout of the remainder of the paper is as follows: In the next section we 
briefly describe the finite lattice method and the nature of the results we have obtained 
therefrom. In section 3 we analyse the data. In section 4 we present a discussion of 
the results. 

2. Series expansions from the finite lattice method 

The definitions and notation follow the usage of I (and II). The standard g-state 
Potts model is defined on a lattice with each site having a 'spin' variable that takes 
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1. Introduction 



This is the third in a series of papers in which we study the critical behaviour of 
the g-state Potts model in both two and three dimensions using series expansions 
derived from the finite lattice method. The first paper (Guttmann and Enting, 
1993a), denoted I hereafter, gave the general expressions used to derive high- and 
low-temperature expansions for the g-state Potts model. In I, series expansions for 
the q = 2 (Ising) case on the simple cubic lattice were analysed. The second paper 
(Briggs et al 1994) denoted II hereafter, presented and analysed series for the bulk 
thermodynamic properties for Potts models on the square lattice for integer q ranging 
from 2 to 10. These were used to develop and test series analysis techniques that 
could distinguish between first-order and continuous transitions. The present paper 
considers the three-state model on the simple cubic lattice. 

After the initial paper by Potts (1952), the model attracted little attention for 
almost two decades. During the 1970's there was greatly renewed interest in the 
model, with new exact results, series studies and renormalisation group calculations 
and applications to phase transitions in surface films. A particular concern at that time 
was the failure of renormalisation group calculations to reproduce the exact results 
for the order of the transition in two dimensions. A review by Wu (1982) described 
much of the work on the Potts model. 

Of even greater interest is the behaviour of the three-dimensional Potts model. As 
noted above, for the q = 2 (Ising) case, the low-temperature series and some high- 
temperature series have recently been extended in I. For 9 = 3, the three-dimensional 
Potts model is of particular interest as it is in the same universality class as the Z(3) 
clock model. This in turn is the centre of SU(3), so it is believed that the effective 
theory for Polyakov loops in finite-temperature d = 4 SU(3) lattice gauge theory should 
be in the same universality class as the three-dimensional three-state Potts model. 
The Polyakov loops are the order parameter for the deconfinement transition that is 
thought to take place in QCD as the temperature is raised (hadronic quarks going over 
to a plasma of free quarks and gluons). The connection then is that if the d = 3, q = 3 
Potts model has a first-order transition then the above SU(3) transition must also be 
first-order. On the other hand if the Potts model has a continuous transition then the 
SU(3) transition could be either first-order or continuous. A number of earlier studies 
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Abstract. The finite lattice method of series expansion has been used to extend low- 
temperature series for the partition function, order parameter and susceptibility of the 
3-state Potts model on the simple cubic lattice to order z 43 and the high-temperature 
expansion of the partition function to order v 21 . We use the numerical data to show 
that the transition is first-order, and estimate the latent heat, the discontinuity in the 
magnetisation, and a number of other critical parameters. 
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